'''
Created on May 19, 2009

@author: mkiyer
'''

from optparse import OptionParser
from veggie.parser.oldelandparser import ElandParser
from veggie.parser.bedparser import get_hpeak_bed_line, get_bed_line
import sys
import re

def convert_elandresult_2_bed(elandfile, outfile=sys.stdout, do_hpeak=False):
    '''
    convert an eland export file to BED file format
    
    elandfile: file handle to eland export file
    outfile: file handle to BED output file
    '''
    # parse file
    for read in ElandParser.parseElandResultFile(elandfile):
        # only process reads with unique best matches
        if read.hasUniqueBestMatch():
            # maps to contaminant file
            if read.chrom == 'newcontam.fa':
                continue
            # maps to ribosomal
            if read.chrom == 'humRibosomal.fa':
                continue
            # maps to mitochondrial
            if read.chrom == 'chrM.fa':
                continue
            # reads mapping to genome
            if (read.chrom.startswith('hs_ref')):        
                # convert strand format (F, R) to BED (+, -)
                strand = '+' if (read.strand == 'F') else '-'
                # convert hs_ref_chrN format to chrN (remove the hs_ref_)
                chrom_name = re.match(r'hs_ref_(\w+).fa', read.chrom).group(1)
                # create BED read format
                if do_hpeak:
                    bedline = get_hpeak_bed_line(chrom_name, 
                                                 str(read.position),
                                                 str(read.position + len(read.seq)),
                                                 strand)
                else:
                    bedline = get_bed_line(chrom_name, 
                                           str(read.position),
                                           str(read.position + len(read.seq)),
                                           strand)
                outfile.write(bedline + '\n')
                continue

def convert_elandexport_2_bed(elandfile, outfile=sys.stdout, do_hpeak=False):
    '''
    convert an eland export file to BED file format
    
    elandfile: file handle to eland export file
    outfile: file handle to BED output file
    '''
    # parse file
    for read in ElandParser.parseElandExportFile(elandfile):
        # only process reads passing filter
        if read.passes_filter:
            # maps to contaminant file
            if read.match_chromosome == 'newcontam.fa':
                continue
            # maps to ribosomal
            if read.match_chromosome == 'humRibosomal.fa':
                continue
            # maps to mitochondrial
            if read.match_chromosome == 'chrM.fa':
                continue
            # TODO: reads matching one of the adaptors?
            # has QC issues
            if read.match_chromosome == 'QC':
                continue
            # non-mapping reads
            if read.match_chromosome == 'NM':
                continue
            # reads with too many multi-maps (e.g. 12:200:255)
            if re.match(r"^[0-9]", read.match_chromosome):
                continue
            # TODO: handle reads with match_chromosome == allMrna.fa
            # reads mapping to genome
            if (read.match_chromosome.startswith('hs_ref')):        
                # convert strand format (F, R) to BED (+, -)
                strand = '+' if (read.match_strand == 'F') else '-'
                # convert hs_ref_chrN format to chrN (remove the hs_ref_)
                chrom_name = re.match(r'hs_ref_(\w+).fa', read.match_chromosome).group(1)
                # create BED read format
                if do_hpeak:
                    bedline = get_hpeak_bed_line(chrom_name, 
                                                 str(read.match_position),
                                                 str(read.match_position + len(read.read_seq)),
                                                 strand)
                else:
                    bedline = get_bed_line(chrom_name, 
                                           str(read.match_position),
                                           str(read.match_position + len(read.read_seq)),
                                           strand)
                outfile.write(bedline + '\n')
                continue
            
            # reads mapping to ucsc genes
            if (read.match_chromosome.startswith('ucsc')):
                gene = read.match_contig
                pos = read.match_position
                strand = '+' if (read.match_strand == 'F') else '-'
                length = len(read.read_seq)
                #print 'ucsc_genes map:', gene, pos, strand, length
                # TODO: for now this util is used for chipseq so ignore
                # the ucsc_genes mappings, but eventually make the eland parser
                # customizable to output different types of information according to
                # user preferences
                continue

if __name__ == '__main__':
    # setup command line parser
    usage = "usage: %prog [options] elandfile"
    optparser = OptionParser()
    optparser.add_option("--hpeak",
                         action="store_true", dest="hpeak", default=False,
                         help="create the pseudo-BED format needed by hpeak")
    optparser.add_option("--format", action="store", dest="format", choices=('export', 'eland_result'), default='export')
    optparser.add_option("-o", "--output",
                         metavar="FILE", help="write output to FILE (default output is to terminal)")
    
    # parse command line
    (options, args) = optparser.parse_args()
    # check options, args
    if len(args) < 1:
        optparser.error('No input file specified')
    
    # get arguments
    elandfilename = args[0]
    
    # setup output
    if options.output:
        outfile = open(options.output, 'w')
    else:
        outfile = sys.stdout
    
    # parse/convert
    if options.format == 'export':
        convert_elandexport_2_bed(open(elandfilename), outfile, options.hpeak)
    elif options.format == 'eland_result':
        convert_elandresult_2_bed(open(elandfilename), outfile, options.hpeak)
